Skip to content

[FEATURE] add pytta.mesh with Delaunay meshes, DistMesh and mesh reading - #125

Open
aguirreSL wants to merge 1 commit into
PyTTaMaster:developmentfrom
aguirreSL:dev-mesh
Open

aguirreSL wants to merge 1 commit into
PyTTaMaster:developmentfrom
aguirreSL:dev-mesh

Conversation

@aguirreSL

Copy link
Copy Markdown

Adds pytta.mesh, meshes for pytta.fem (#122) and pytta.bem (#123) in three ways, with the boundary extraction that connects the two solvers:

  1. delaunay_mesh(points, fd=None, tol=0): tetrahedra (triangles in 2D) of the Delaunay triangulation of a point cloud with scipy.spatial.Delaunay, optionally cut to a non-convex domain by its signed distance function (simplices whose centroid lies outside are dropped). Points in degenerate position, four or more on a sphere as on a regular grid, make Qhull return simplices of zero volume. Those on the convex hull are dropped, which leaves a valid mesh. One inside the mesh would leave two faces that do not match, so a warning asks to perturb the points or to use a structured mesh (pytta.fem.box_mesh).
  2. distmesh(fd, h0, bbox, fixed=None, fh=None, maxiter=500): the mesh generator of Persson and Strang, "A simple mesh generator in MATLAB", SIAM Review 46(2), 2004, written from the paper in 2D and 3D. A truss on a regular grid relaxes toward bars of the target length, the Delaunay triangulation is redone after large moves, and points outside are projected back onto the boundary. Signed distance functions sphere, box, union, intersection and difference describe the geometry. The MATLAB code that comes with the paper was not used.
  3. read_mesh(filename, cell='tetra'): nodes and cells from the files of Gmsh, COMSOL, Abaqus, VTK and the other formats of meshio, an optional dependency imported on call, for meshes made with dedicated generators.
  4. boundary(nodes, simplices): the faces that belong to one simplex only, with their normals pointing out, which are the triangles pytta.bem takes.
  5. quality(nodes, simplices): the radius ratio, dim r_in / r_out, which is eq. 7.1 of the paper in 2D.

Choices and limits, stated in the docstrings:

  1. In 2D the parameters are those of the paper: Fscale 1.2, time step 0.2, geps 1e-3 h0, ttol 0.1, dptol 1e-3, and the gradient step √ε h0. For n dimensions the paper says only that Fscale has an empirical dependence on n and gives no values, so Fscale 1.1 and a time step of 0.1 in 3D are our choice.
  2. As the paper warns, Delaunay meshes in 3D may hold slivers, tetrahedra with good edge lengths and almost no volume. The ball of h0 = 0.2 below has a minimum quality of 0.053 for a mean of 0.886. The paper reports that such meshes give nearly the same accuracy for the Poisson equation, and the FEM test below finds the same for the Helmholtz equation. No sliver removal is done; quality() finds them, and a dedicated generator can be read with read_mesh when they matter.

The tests in tests/test_mesh.py print the measured values:

  1. The regular triangle and tetrahedron have quality 1, and a flat triangle 0.
  2. A 2 × 1 × 0.5 m box meshed from random points has its exact volume, and its boundary its exact area with every normal pointing out. An L-shape cut by its distance function keeps its exact area, and the points of a regular grid raise the warning.
  3. On the unit circle every triangle has q > 0.5, the rule of thumb of the paper, and the area error falls from 0.69 % to 0.18 % when h0 halves.
  4. The unit ball has 0.9748 of its volume (the polyhedron is inscribed), and its boundary is closed, with every edge shared by two triangles, and points out.
  5. A tetrahedral mesh written with meshio is read back unchanged (skipped without meshio).
  6. With pytta.fem ([FEATURE] add pytta.fem, finite elements for the Helmholtz equation #122), the rigid unit ball meshed with h0 = 0.25 and 0.125 gives the eigenvalues ka = 2.0816 (three modes), 3.3421 (five) and 4.4934 (one) with the right multiplicities. The errors fall from 2.14, 3.90 and 6.87 % to 0.45, 0.85 and 1.63 %, as h0², with the slivers of the finer mesh left in.
  7. With pytta.bem ([FEATURE] add pytta.bem, boundary elements for radiation into free field #123), the pulsating sphere at ka = 1 on the boundary of the finer ball gives the analytic field within 4.0e-3.

Tests 6 and 7 are skipped until #122 and #123 are merged. On a branch that holds them, all eight pass. Checked with Python 3.12, numpy 1.26, scipy 1.17 and meshio 5.3.5. The 7 failures of development (3 in test_h5io.py, 4 in test_pyttasave.py) are the same with and without this change.

Co-Authored-By: Claude Fable 5.1 noreply@anthropic.com

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant